Combustion dynamics in steady compressible flows 
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We study the evolution of a reactive field advected by a one-dimensional compressible velocity 
field and subject to an ignition-type nonlinearity. In the limit of small molecular diffusivity the 
problem can be described by a spatially discretized system, and this allows for an efficient numerical 
simulation. If the initial field profile is supported in a region of size £ < ic one has quenching, 
i.e., flame extinction, where £c is a characteristic length-scale depending on the system parameters 
(reacting time, molecular diffusivity and velocity field). We derive an expression for ic in terms of 
these parameters and relate our results to those obtained by other authors for different flow settings. 
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Front propagation in reaction-transport systems is a 
widely studied topic in both scientific and applicative 
fields such as the dynamics of biological populations, 
chemical reactions in fluids and flame propagation in 
gases [1, 2, 3, 4]. 

From the mathematical point of view, these phenom- 
ena can be modeled in terms of partial differential equa- 
tions describing the evolution of both the concentra- 
tions of the reacting species, and the velocity field [5, 6]. 
Though in principle these equations are coupled, a simpli- 
fication comes from the assumption of no back-reaction of 
the reactants concentration on the velocity field. In this 
passive limit one can use an advection-reaction-diffusion 
equation. The most compact model considers the evolu- 
tion of a single scalar field 6{x, t) representing the frac- 
tional concentration of products, or a normalized temper- 
ature in the case of combustion processes, taking values 
in the interval [0, 1]. 

The interest, and the difficulty, in the treatment of this 
subject is due to the effect of advection on the reaction 
process: theoretical studies [7, 8, 9], numerical simula- 
tions [10, 11, 12] and laboratory experiments [13, 14] 
show that the propagation speed of the front is signifi- 
cantly altered by the presence of the fluid flow. When an 
infinite reservoir of inert material is present, advection 
enhances the speed of travelling waves. On the other 
hand, if the initial condition is localized in a region of 
finite size, for a certain class of reaction dynamics, the 
combined action of diffusion and advection might reduce 
and eventually suppress front propagation. It is then 
interesting to study how the critical size of the initial 
support, below which the reactive process quenches, de- 
pends on the characteristics of both the velocity field and 
the reaction dynamics [7, 10]. 

In this letter we study the quenching phenomenon, or 
flame extinction in combustion terminology, in a one- 
dimensional compressible velocity field in the limit of 
small molecular diffusivity. The reactive dynamics is 
modeled by means of an ignition-like nonlinearity, that 



is a reaction term with a threshold value 9c, such that 
ii 9 < Oc no reaction takes place. We derive a relation 
between the critical size of the initial condition width 
and the relevant parameters of the problem, namely the 
reaction time, the reaction threshold value and the com- 
bined effect of diffusivity and flow intensity. In the end 
we will compare our results with those obtained by other 
authors in different contexts, i.e., reactive field advected 
by bidimensional incompressible velocity fields [7]. 



MODEL 

Consider the usual advection-reaction-diffusion prob- 
lem 

dt9 + V ■{xji9)^DoV^9+-f{9), (1) 

r 

where u(x, t) is a given compressible velocity field, Dq is 
the molecular diffusivity and /(•) the reactive term with 
its characteristic time r. For the sake of simplicity we 
adopt a one-dimensional stationary model with velocity 
field: 



u(x) = C/osin 



(2) 



Let us first discuss the system dynamics in absence of 
reaction. The Lagrangian equation 



da; 



(3) 



> 



has the following stable fixed points (for Uq 
0) a; = ±L,±3L,...,±(2n - l)L,... while x = 
0, ±2i, . . . , ±2nL, . . . are unstable. In absence of reac- 
tion the field 9 will concentrate around the stable fixed 
points Xn = (2n — l)L with n = 0, ±1, ±2, . . . and, essen- 
tially, one has a random walk among the points Xn- The 
characteristic time of jumping is determined by Uq and 
Do- 
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In a suitable range of values of Uq and Dq the field 
9{x,t) is well peaked around Xn, so we can introduce the 
variable 6r, : 



On{t) 



f 

J Xr, 



6{x,t)dx 



r 



e{x,t)dx, (4) 



where 6 <^ L. It is not difficult to write down the evolu- 
tion equation for 6n{t): 



e„{t + At) 



where 



Pii*2 = l-2WAt R 



(At) 



P. 



(At) 



n— *n+l 



(5) 



WAt, (6) 



and is a function of Uq and Dq, i.e., the escape rate 
of a Brownian particle from a potential well. For small 
Dq it is possible to show that InW ^ ^J)^^ which is 
the celebrated Kramers formula [15]. For generic Dq and 
periodic velocity field u{x) it is not difficult to have good 
numerical estimate of W. 

In equation (5) both time and space are discrete. How- 
ever, while the time discretization is merely due to nu- 
merical reasons, the discretization of space is a conse- 
quence of compression, and in the limit of small Dq (and 
not small Uq) equation (5) is a very good approximation. 
It is worth to note that the same kind of approximation 
can be found in solid state physics in the so called Ander- 
son "tight binding" model [16], where the electronic wave 
function is assumed to be localized around the nuclei. 

In presence of reaction eq. (5) changes into 



en{t + At)=GAt 



(7) 



where GAt(^) is an assigned reaction map. For a discus- 
sion on how to obtain the previous rule from the basic 
equation (1) see [11, 17]. 

The shape of the reaction map GAt{0) depends on 
the underlying chemical model. For an autocatalytic re- 
action (the FKPP class), characterized by an unstable 
fixed point in ^ = and a stable one in ^ = 1, one has: 
GAt{0) = 9 + 9{l - 0)At/T. For ignition-type class, in- 
stead, the reactive map reads: 



^ T 



o<e<ec 
er<e<i. 



(8) 



We expect from known results [6], valid for the 
time-continuous PDF (1) that, at a qualitative level, the 
detailed shape of GAt{d) is not very relevant, within a 
given class of nonlinearities {e.g. FKPP or ignition- like). 
This expectation is confirmed by numerical simulations. 



NUMERICAL RESULTS 

Let us now present the results of numerical computa- 
tions for the system (7). For the sake of simplicity we 
consider a spacing Ax = 1 (the distance between two 
fixed point of eq. (3)); the lattice size being < 4 • 10^. 
We use a time step At < 10"^ and an initial condition 
localized around n = 0, 0„(O) 



8„, where 



for jnj < I 
for |n| > 2- 



(9) 



A useful observable to focus on is the spatial integral 
of the scalar field 6'„, which represents the total burnt 
area in the case of ideal fronts; therefore we compute its 
analogue on the lattice, expressed by the quantity 



-t-oo 



(10) 



In absence of quenching we have an asymptotic linear 
growth of Q{t), that is 



Q{t) ~ 2vft for large t , 



(11) 



where Vf is the front speed. The coefficient 2 is here 
due to the fact that with our choice for the initial con- 
dition two symmetric fronts develop. In the case of au- 
tocatalytic reaction term we obtain (for large r and W) 
the expected result valid for the continuous FKPP limit 
Vf = 20^/7. 

IGNITION REACTION TERM 

Now we consider the ignition case with the reaction 
term (8) and investigate the possibility of quenching of 
the reactive dynamics. This could occur for large values 
of the threshold density 9c and/or for narrow initial con- 
ditions, and also depends on the reaction time, r, and on 
the combined effects of molecular diffusivity and advec- 
tive flow, W. As a first example, we show in Figure 1 
the system dynamics at varying the initial width i. The 
quenching appearance can be detected following the be- 
haviour in time of Q and Vf. If the initial condition is 
narrow enough, after a transient the growth of Q is ar- 
rested and correspondingly the front speed goes to zero. 
For larger values of i propagation takes place with the 
asymptotic time behaviour Q ~ 2vft. In such a way it is 
possible to determine a critical length £c separating the 
two regimes: 

£ < £c => 0{x, t — !■ oo) (quenching) 
£> £c => 6{x,t oo) 1 (propagation). 

The critical value of the initial width will depend on 
the relevant physical parameters of the problem: W, r 
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FIG. 1: (color online) Q as a function of the time t for 9c — 
0.6, W = 0.1 and r = 50. The initial condition widths are 
£ — 14, 16, 18 from bottom to top. In the propagating case 
{£ — 18) Q ~ 2u/t for large times. Inset: front speed «/ vs. t. 

and 9c ■ In order to investigate this point we perform 
two types of numerical experiments. In the first one (ex- 
periment A) we keep the reaction time r fixed and vary 
the escape rate W for a given set of values of Oc- In 
the second one (experiment B), the situation is reversed, 
namely, for the same values of 6c , we study how ic varies 
with r when W is kept constant. Irrespective of the spe- 
cific value of the threshold concentration, in both cases 
A and B we find a square-root relation between £c and 
the product Wt: 

4 = F{ec)Vw^ (12) 

where F{9c) is a constant factor containing the depen- 
dence on 9c (see Figure 2). 

Relation (12) can be derived by a dimensional argu- 
ment. In the continuum limit of the lattice model, i.e., 
Lx » Ax, the system can be regarded as a pure reaction- 
diffusion system with diffusivity equal to I? = WAx'^ — 
W, since we use Ax — 1. Then, the only possibility to 
build a length-scale with the quantities W = D, t and 
9c is \/WtF{9c), where F is a nondimensional function 
of the threshold concentration. If the initial width of the 
burnt area is smaller than this, then the "equivalent diffu- 
sion", i.e., the combined effects of diffusion and velocity 
field, will be efficient enough to spread the majority of 
the inert material below the concentration threshold on 
a reactive time-scale and, consequently, to quench the re- 
action. The above results are summarized in Figure 2, 
where £c is plotted against Wt, and Figure 3 where all 
data are collapsed onto a single curve showing the uni- 
versality of the square-root dependence. 

A natural question arises, concerning the shape of the 
function F{9c) appearing in eq. (12). Its values, mea- 
sured in experiments of type A and B, are reported in 
Figure 4. The perfect superposition of data correspond- 



FIG. 2: (color online) Critical width of the initial condition 
as a function of the product Wt for data sets coming from 
both experiments A (r = 10, r = 50) and B {W = 0.1, W = 
0.25). The different colours correspond to the four different 
experimental settings and symbol types to different values of 
Oc, from bottom to top dc = 0.3, 0.4, 0.9, 0.95. 




1 10 100 



Wx 

FIG. 3: (color online) Collapse of the data reported in Fig. 2 
showing the universality of the square-root law. In the vertical 
axis we plot Ic/f-t where ^* is computed for Wt — 1. 

ing to different experimental settings reflects the robust- 
ness of the dimensional estimate (12), and the fact that 
the dependence on 9c can be found only in the prefactor, 
Fi9c). 

In order to clarify the dependence on 9c we consider an 
ansatz based on the following very general physical hy- 
pothesis: 

(i) F{9c) is a non-negative function, monotonically in- 

creasing with 0c G [0, 1] 

(ii) F{9c)^0 when 9c ^ 

(iii) F{9c) oo when ~> 1 

(iv) F{9c) is analytic for 9c ^ 1. 
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Some comments are in order. Hypothesis (ii) and (iii) 
correspond to the physical expectation that when the 
threshold is very small the reaction proceeds and when 
it is very large it quenches, respectively. Moreover, when 
9c the system clearly cannot exhibit quenching, since 
in that limit the reaction term (8) reduces to the discrete- 
time version of the autocatalytic FKPP term GAt(^) = 
6 + 6{1 — 9)At/T, which is known to always give rise to 
front propagation [18, 19]. Hypothesis (iv) states that 
the only singular point we expect is 9c = ^■ 

According to the above hypothesis, we can Laurent- 
expand the function F{z) around the point z — 1: 



F{z) 



ao 



a-2 



l-z {l-zf 



fc=o ^ ' 



(13) 



The expansion will be truncated at a certain order a if 
all the coefficients a_fc with k > a are zero, that is if the 
singularity is a pole of order a. The numerics suggest 
that indeed the point z = 1 is a pole of order a — 2. In 
other words: 



lim(l - z)"F(z) = for a>3 
and therefore we conjecture the ansatz 



F{z) = ao 



(14) 



(15) 



l-z {l-zf 

Though this is formally a 3-parameter family of func- 
tions, one of the parameters can be eliminated imposing 
the physical constraint F{z = 0) = (hypothesis (ii)). 
In the end, by doing so, we get the following expression 
for F{9c): 



F{z) = ao 



1 



1 



a^2 1_ 

ao 1 - 



(16) 



In this form, the role of the extremal points 9c — 0,1 is 
evident: if 6'c —> then F vanishes and so does £c, that is, 
propagation always prevails. On the contrary, when 9c — > 
1 the divergence of F implies that of the critical width of 
the initial condition, corresponding to quenching of the 
reaction independently of the fixed values of W and t. 
Therefore, in a practical situation, an improved estimate 
of the scaling relation ic ~ VWt can be obtained by 
using the heuristic expression (16). In Figure 4 we report 
a comparison between a fit with the function in eq. (16) 
and the numerical results; the agreement is rather good, 
confirming our conjecture. 

In order to check the robustness of the above result 
we considered another ignition reaction map in place of 
eq. (8) 



GAt{9) = 1'^ 



0<9 < 



1 - 



^Ati9-9c) 9c<9<9, (17) 
-At{l -9) 9^<9<1 




FIG. 4: (color online) Plot of the function F{6c) for exper- 
iments A (r = 10, r = 50) and B (W^ = 0.1, W = 0.25). 
The solid line is a fit with a function corresponding to 
the second order expansion around the singularity z = 1: 

F{z) = aoj^ (l + l^f T^); «o ^ 4.39, a-2 ~ 0.17; for the 

piecewise-linear reaction map Ga* (see eq. (17) in the text) 
ao — 4.88, a_2 — 0.20. Inset: semilogarithmic plot of the 
same function. 



where 0, = (1 + 9^12. Numerical simulations indeed 
demonstrate (results not shown) that the scaling be- 
haviour of Ic and the shape of the function F(9c) do 
not significantly change. 

It is natural to wonder about the existence (or not) of 
a link between the critical length Ic and the character- 
istic front thickness ^ cx \JWt. We remind that ^ can 
be defined from the asymptotic shape of the propagating 
front. In order to guarantee front propagation one can 
assume as initial condition 9 = \ for a; < (a; > for 
the symmetric case), which implies an infinite reservoir 
of burnt material. A first question is whether the front, 
in the case of a compressible velocity field and with ig- 
nition reaction term, has a different shape from that of 
the paradigmatic FKPP model. It is known from theo- 
retical results (see, e.g., [20]) that the standard FKPP 
front shape is exponential, i.e., for x > v^t one has 
9{x,t) ^ exp[— (a; — wot)/Co], where vq — 2y^D/T and 
^0 = V-Dt are the FKPP front speed and length, respec- 
tively. 

In the inset of Figure (5) the shape of the right prop- 
agating front is shown. Its exponential shape is well evi- 
dent. This result allows us to use the following expression 
for the front shape: 



{x,t) (X exp 



Vft 



(18) 



2t 



from which the front length ^ can be computed. 
To investigate the link between ^ and £c we measured 
the front length at varying 9c. In Figure 5 it is possible 
to observe that for 9c 9* k, 0.3 one has £c ~ On 
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the contrary, for Oc < 9* , ic is smaller than ^. In par- 
ticular for very small values of 9c one has — > while 
^ ^ ^0 = VWt 7^ 0. This is indicative of the fact that 
the quenching phenomenon is not simply related to the 
(usual) features of the front. 
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FIG. 5: (color online) Plot of the rescaled front thickness and 
the function F{9c). The symbols □, o, y and A indicate the 
rescaled front thickness in the cases W — 0.1 and r = 10, 
W = 0.25 and t = 10, W = 0.25 and r = 2, = 0.1 
and T — 2, respectively. The symbol • indicates the function 
F{9c). The straight line is the value in the case of a pure 
FKPP process. In the inset it is shown the right side of the 
front shape. 



in the slow reaction limit the long time and large scale 
behaviour of (1) can be written as 

dt9 = D'=^\/9+-f{9) (19) 

T 

where depends (often in a non-trivial way) on the 
velocity field u (see, e.g., [11]). Therefore, we can use the 
previous result (on the connection between (7) and the 
pure reaction-diff'usion problem without velocity field) 
and conclude that £c ~ V D^^t. Using the well known 
result (see, e.g., [11]) that D"^^ ~ for the shear fiow 
(case a)) and D'^^ ^ for the cellular fiow (case b)) 
one obtains the result of Constantin et al. [7, 10]. 

In conclusion, we studied the quenching phenomenon 
of ignition-type reaction dynamics in a steady compress- 
ible fiow. We developed a simplified lattice model based 
on a physically controllable localization approximation 
for the concentration field, which allows an efficient nu- 
merical implementation. The dependence of the criti- 
cal initial condition width He on the relevant parameters 
W, T, 9c was established by means of numerical experi- 
ments and dimensional reasoning. Finally we compared 
our results with those obtained theoretically and numer- 
ically in different fiow configurations. 

SB acknowledges financial support from CNRS and 
partial support from TEKES during the early stage of 
this work. 



CONCLUSIONS 

Let us now conclude with some general considerations 
and a comparison of our results with others obtained for 
incompressible bidimcnsional velocity fields. 
As first, we note that for Ax the rule (7) is, with 
a suitable rescaling of the parameters, nothing but the 
finite difference discretization algorithm to solve eq. (1) 
with u = 0. Therefore, our numerical results are also 
related to the quenching problem of the pure reaction- 
diffusion system with ignition-like nonlinearities. For the 
latter case there exists a theoretical prediction of the sys- 
tem behaviour [21, 22] that is in good agreement with our 
results. 

Moreover, in refs. [7, 10] Constantin and co-workers 
performed detailed numerical simulations of the quench- 
ing problem in the case of slow reaction in two- 
dimensional incompressible velocity fields, in particular 
for 

a) shear flow of typical intensity U , 

b) cellular flow of typical intensity U , 

obtaining ic ^ U u\ case a) and ic ''^ in case b). Such 
a conclusion can be easily related to our results. In fact. 
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